slope <- map_dfr(rigal_trends, ~.x %>%
  select(siteid, linear_slope),
.id = "response"
)
slope_df <- slope %>%
  pivot_wider(names_from = "response", values_from = "linear_slope")

4 Comparison rigal, lm

rigal_trends_df <-  map2_dfr(
  rigal_trends, names(rigal_trends),
  ~.x %>% mutate(variable = .y)
  ) %>%
  select(variable, siteid, linear_slope) %>%
  pivot_wider(names_from = "variable", values_from = "linear_slope") %>%
  mutate(type = "rigal")
slope_comp <- list()
slope_comp$rigal <- rigal_trends_df[order(rigal_trends_df$siteid), ]
slope_comp$simple_lm <- slope_df[order(slope_df$siteid), ]
stopifnot(all(slope_comp$rigal$siteid == slope_comp$simple_lm$siteid))
slope_comp$diff <- tibble(
  siteid = slope_comp$rigal$siteid)
for (i in var_temporal_trends) {
  slope_comp$diff[[i]] <- slope_comp$rigal[[i]] - slope_comp$simple_lm[[i]]
}
old_par <- par()
par(mfrow = c(3, 2))
for (i in var_temporal_trends) {
  plot(
    slope_comp$rigal[[i]] ~
      slope_comp$simple_lm[[i]],
    ylab = "Rigal (polynomial degree 2)",
    xlab = "Simple LM"
  )
  abline(0, 1)
  title(i)
}

par(old_par)
#> Warning in par(old_par): le paramètre graphique "cin" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "cra" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "csi" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "cxy" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "din" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "page" ne peut être changé

# Compare Rigal and simple lm for few sites
plot_rigal_raw_data <- function(
  site = NULL,
  response = NULL,
  dataset = NULL,
  rigal_trends = NULL
  ) {

  raw_data <- dataset %>%
    filter(siteid == site)
  rigal_resp <- rigal_trends[[response]]
  rigal_resp_site <- rigal_resp[rigal_resp$siteid == site, ]
  plot(raw_data[[response]]~raw_data[["year"]], ylab = response, xlab = "Year")
  abline(rigal_resp_site$intercept, rigal_resp_site$linear_slope, col
  = "red")
  abline(lm(raw_data[[response]]~raw_data[["year"]]), col = "green")
  title(main = site)
}
plot_rigal_raw_data(site = "S2672", response = "log_total_abundance", dataset =
analysis_dataset, rigal_trends = rigal_trends)

slope_comp$diff <- slope_comp$diff %>%
  arrange(desc(abs(log_total_abundance)))
old_par <- par()
par(mfrow = c(3, 2))
for (i in slope_comp$diff[1:6,]$siteid) {
  plot_rigal_raw_data(
    site = i,
    response = "log_total_abundance",
    dataset = analysis_dataset,
    rigal_trends = rigal_trends
  )
}

par(old_par)
#> Warning in par(old_par): le paramètre graphique "cin" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "cra" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "csi" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "cxy" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "din" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "page" ne peut être changé

5 Polynomial

classification_station <- map2(rigal_trends, names(rigal_trends),
  ~.x %>%
  ggplot(aes(x = shape_class)) +
  geom_bar() +
  labs(x = "Trajectory classification", y = "Nb station",
    title = get_var_replacement()[.y]) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
)


plot_grid(plotlist = classification_station)

p_classif_abun_rich <- plot_grid(
  plotlist = classification_station[
    names(classification_station)[str_detect(names(classification_station),
  "abundance|richness|species")]])
save_plot(
  here("doc", "fig", "classif_trends_abun_rich.png"),
  p_classif_abun_rich
)
p_classif_hill_evenness <- plot_grid(
  plotlist = classification_station[
    names(classification_station)[str_detect(names(classification_station),
  "shannon|simpson|evennes")]])
save_plot(
  here("doc", "fig", "classif_trends_hill_evenness.png"),
  p_classif_hill_evenness
)
p_classif_turnover <- plot_grid(
  plotlist = classification_station[
    names(classification_station)[names(classification_station) %in%
      c("jaccard", "horn", "chao", "hillebrand", "total", "appearance",
        "disappearance")]])
save_plot(
  here("doc", "fig", "classif_trends_classif_hill_evenness.png"),
  p_classif_turnover
)
ti <- map_dfr(rigal_trends, ~tabyl_df(x = .x, group = "direction"),
  .id = "response"
)
ti %>%
  filter(direction != "Total") %>%
  select(response, direction, percent) %>%
  pivot_wider(names_from = "direction", values_from = "percent")
#> # A tibble: 21 × 4
#>    response            stable increase decrease
#>    <chr>               <chr>  <chr>    <chr>   
#>  1 total_abundance     75.8%  14.3%    9.9%    
#>  2 log_total_abundance 73.3%  16.5%    10.2%   
#>  3 species_nb          77.7%  13.2%    7.4%    
#>  4 log_species_nb      76.2%  13.2%    6.6%    
#>  5 chao_richness       79.8%  10.6%    8.2%    
#>  6 chao_shannon        76.5%  12.1%    7.7%    
#>  7 chao_simpson        78.4%  10.5%    7.4%    
#>  8 chao_evenness       82.8%  7.5%     5.9%    
#>  9 jaccard             76.0%  3.5%     20.4%   
#> 10 horn                75.9%  4.4%     17.6%   
#> # … with 11 more rows
#3.3   Log-linear model:logYi=α+βXi+εiIn

#the log-linear model, the literal interpretation of the estimated
#coefficientˆβis that a one-unitincrease inXwill produce an expected increase in
#logYofˆβunits. 
#In terms ofYitself, this
#meansthat  the  expected  value  ofYis  multiplied  byeˆβ.   So  in  terms  of
#effects  of  changes  inXonY(unlogged):•Each 1-unit increase inXmultiplies the
#expected value ofYbyeˆβ.•To compute the effects onYof another change inXthan an
#increase of one unit, call thischangec, we need to includecin the exponent.  The
#effect of ac-unit increase inXis tomultiply the expected value ofYbyecˆβ. So the
#effect for a 5-unit increase inXwould bee5ˆβ.•For small values ofˆβ,
#approximatelyeˆβ≈1+ˆβ. We can use this for the following approxima-tion for a
#quick interpretation of the coefficients:  100·ˆβis the expected percentage
#changeinYfor a unit increase inX.  For instance forˆβ=.06,e.06≈1.06, so a 1-unit
#change inXcorresponds to (approximately) an expected increase inYof 6%
#https://data.library.virginia.edu/interpreting-log-transformations-in-a-linear-model/
hist_linear_plot <- map2(rigal_trends, names(rigal_trends),
  ~.x %>%
  ggplot(aes(x = linear_slope)) +
  geom_histogram() +
  labs(x = "Linear slope", y = "Number of station",
    title = get_var_replacement()[.y]) +
  ylim(c(0, 1000))
)
plot_grid(plotlist = hist_linear_plot)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).

type_response <- list(
  abun_rich = str_detect(names(classification_station),
    "abundance|richness|species"),
  hill_evenness = str_detect(names(classification_station),
    "shannon|simpson|evennes"),
  turnover = names(classification_station) %in%
    c("jaccard", "horn", "chao", "hillebrand", "total", "appearance",
      "disappearance")
)

map2(type_response, names(type_response),
  function(x, y) {
  p <- plot_grid(
    plotlist = hist_linear_plot[
      names(hist_linear_plot)[x]
      ]
  )
  save_plot(
    here("doc", "fig", paste0("hist_linear_slope_", y, ".png")),
    p
  )

  }
)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> $abun_rich
#> NULL
#> 
#> $hill_evenness
#> NULL
#> 
#> $turnover
#> NULL

5.2 Sensibility analysis

tar_load(c(rigal_trends_avg3y))
names(rigal_trends_avg3y) <- var_temporal_trends
tar_load(c(analysis_dataset_avg3y, measurement_avg3y, measurement))
tar_load(c(turnover_avg3y_c))
tar_load(c(vegdist_turnover_c))
tar_load(c(hillnb_avg3y))
slope_avg3y <- map_dfr(rigal_trends_avg3y[
  ! names(rigal_trends_avg3y) %in% c("jaccard_dis", "nestedness", "turnover")], ~.x %>%
  select(siteid, linear_slope),
.id = "response"
)
  • Number of sites that are in the different datasets:
u_site <- map(list(classic = slope, avg_3y = slope_avg3y), ~unique(.x$siteid))
sum(u_site[[1]] %in% u_site[[2]])
#> [1] 2678
sum(u_site[[2]] %in% u_site[[1]])
#> [1] 2678
comparison_baseline <- rbind(
  slope %>%
    mutate(type = "classic") %>%
    filter(siteid %in% u_site[["avg_3y"]]),
  slope_avg3y %>%
    mutate(type = "avg_3y")
  ) %>%
pivot_wider(names_from = type, values_from = linear_slope)
  • It is not too bad, but we can see some differences:
    • slope negatively biaised for log_species_nb in 3 years baseline
    • Bunch of trends in appearance/disappearance trends disappears with
      3 years baseline
    • But note that in three years baseline, there are three points less data
comparison_baseline %>%
  ggplot(aes(x = classic, y = avg_3y)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1) +
  facet_wrap(~response, scale = "free")
#> Warning: Removed 8034 rows containing missing values (geom_point).

5.3 Reproducibility

Reproducibility receipt

## datetime
Sys.time()
#> [1] "2022-01-12 18:11:45 CST"

## repository
if(requireNamespace('git2r', quietly = TRUE)) {
  git2r::repository()
} else {
  c(
    system2("git", args = c("log", "--name-status", "-1"), stdout = TRUE),
    system2("git", args = c("remote", "-v"), stdout = TRUE)
  )
}
#> Local:    main /home/alain/Documents/post-these/isu/RivFishTimeBiodiversityFacets
#> Head:     [0145794] 2022-01-11: Fix shitty R default: drop = FALSE for matrix of one column

## session info
sessionInfo()
#> R version 4.0.5 (2021-03-31)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Debian GNU/Linux 10 (buster)
#> 
#> Matrix products: default
#> BLAS:   /home/alain/.Renv/versions/4.0.5/lib/R/lib/libRblas.so
#> LAPACK: /home/alain/.Renv/versions/4.0.5/lib/R/lib/libRlapack.so
#> 
#> locale:
#>  [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=fr_FR.UTF-8        LC_COLLATE=fr_FR.UTF-8    
#>  [5] LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=fr_FR.UTF-8   
#>  [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] colorspace_2.0-0        future_1.21.0           vegan_2.5-7            
#>  [4] lattice_0.20-41         permute_0.9-5           codyn_2.0.5            
#>  [7] janitor_2.1.0           cowplot_1.1.1           rnaturalearthdata_0.1.0
#> [10] rnaturalearth_0.1.0     sf_1.0-4                rmarkdown_2.11         
#> [13] scales_1.1.1            kableExtra_1.3.1        here_1.0.1             
#> [16] lubridate_1.7.9.2       magrittr_2.0.1          forcats_0.5.1          
#> [19] stringr_1.4.0           dplyr_1.0.4             purrr_0.3.4            
#> [22] readr_2.1.1             tidyr_1.1.2             tibble_3.1.6           
#> [25] ggplot2_3.3.3           tidyverse_1.3.0         tarchetypes_0.3.2      
#> [28] targets_0.8.1           conflicted_1.1.0       
#> 
#> loaded via a namespace (and not attached):
#>  [1] ellipsis_0.3.2     class_7.3-18       rgdal_1.5-28       rprojroot_2.0.2   
#>  [5] snakecase_0.11.0   fs_1.5.1           rstudioapi_0.13    farver_2.0.3      
#>  [9] listenv_0.8.0      fansi_0.5.0        xml2_1.3.2         codetools_0.2-18  
#> [13] splines_4.0.5      cachem_1.0.4       knitr_1.36         jsonlite_1.7.2    
#> [17] broom_0.7.4        cluster_2.1.1      dbplyr_2.1.0       rgeos_0.5-5       
#> [21] compiler_4.0.5     httr_1.4.2         backports_1.2.1    assertthat_0.2.1  
#> [25] Matrix_1.3-2       fastmap_1.1.0      cli_3.1.0          s2_1.0.7          
#> [29] htmltools_0.5.1.1  tools_4.0.5        igraph_1.2.6       gtable_0.3.0      
#> [33] glue_1.5.1         wk_0.5.0           Rcpp_1.0.6         cellranger_1.1.0  
#> [37] jquerylib_0.1.3    vctrs_0.3.8        nlme_3.1-152       xfun_0.28         
#> [41] globals_0.14.0     ps_1.6.0           rvest_0.3.6        lifecycle_1.0.1   
#> [45] MASS_7.3-53.1      hms_1.1.1          parallel_4.0.5     yaml_2.2.1        
#> [49] memoise_2.0.0      sass_0.3.1         stringi_1.7.6      highr_0.9         
#> [53] corrplot_0.84      e1071_1.7-4        rlang_0.4.12       pkgconfig_2.0.3   
#> [57] evaluate_0.14      labeling_0.4.2     processx_3.5.2     tidyselect_1.1.1  
#> [61] parallelly_1.23.0  bookdown_0.24      R6_2.5.1           generics_0.1.0    
#> [65] DBI_1.1.1          pillar_1.6.4       haven_2.3.1        withr_2.4.3       
#> [69] mgcv_1.8-34        units_0.6-7        sp_1.4-5           modelr_0.1.8      
#> [73] crayon_1.4.2       KernSmooth_2.23-18 utf8_1.2.2         tzdb_0.2.0        
#> [77] grid_4.0.5         readxl_1.3.1       data.table_1.13.6  git2r_0.29.0      
#> [81] callr_3.7.0        reprex_1.0.0       digest_0.6.27      classInt_0.4-3    
#> [85] webshot_0.5.2      munsell_0.5.0      viridisLite_0.3.0  bslib_0.2.4